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EXECUTIVE  SUMMARY 

The  objective  of  the  project  was  to  develop  a  multiscale  computational  model  capable  of 
predicting  the  evolution  of  matrix  cracking,  delamination,  and  fiber  cracking  in  viscoelastic 
composite  structures  subjected  to  ballistic  impact.  The  model  is  three  dimensional  and 
computational  in  nature,  utilizing  the  finite  element  method,  and  this  model  is  being 
implemented  to  the  explicit  code  DYNA3D.  Crack  growth  is  simulated  via  the  cohesive  zone 
model  currently  under  development  by  the  author.  The  cohesive  zone  model  for  predicting 
damage  evolution  in  laminated  composite  plates  is  cast  within  a  three  dimensional  continuum 
finite  element  algorithm  capable  of  simulating  the  evolution  of  matrix,  fiber,  and  delamination 
cracking  in  composite  structures  subjected  to  ballistic  impact.  Cracking  on  vastly  differing 
length  scales  is  accounted  for  by  employing  global-local  techniques,  with  appropriate  damage 
dependent  homogenization  techniques  introduced  to  bridge  the  disparate  scales.  Finally, 
simplified  generic  example  problems  were  solved  analytically  and  compared  to  computational 
results  obtained  with  the  model  as  a  means  of  model  verification.  A  damage  constitutive  model 
for  polymer-matrix  composite  materials  was  developed  and  implemented  into  a  commercially 
available  finite  element  package. 


TECHNICAL  DISCUSSION 


Preliminary  Comments 

While  accurate  models  have  been  developed  for  predicting  impact  response  of  metals  and 
ceramics,  which  are  normally  crystalline  in  nature,  models  have  not  been  developed  for 
polymeric  materials,  which  have  the  potential  to  decrease  both  mass  and  cost  of  such  strategic 
defense  weapons.  This  is  due  to  the  different  physical  nature  of  fracture  in  polymeric  media,  as 
well  as  the  multiple  length  scales  on  which  this  fracture  occurs.  The  first  of  these  is  addressed 
by  a  micromechanically  based  viscoelastic  cohesive  zone  model  developed  by  the  principal 
investigator  and  coworkers,  and  the  second  is  solved  via  the  multiscale  damage  dependent 
algorithm  described  herein.  These  two  aspects  of  the  research  comprise  a  unique  approach  to 
this  problem.  The  cohesive  zone  model  is  unlike  any  others  in  existence,  in  that  it  requires  only 
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well  defined  experiments  at  the  microscale.  The  multiscale  algorithm  utilizes  damage  dependent 
homogenization  principles  for  viscoelastic  media  that  have  been  recently  developed  by  the 
author  and  coworkers,  an  aspect  which  has  not  been  correctly  incorporated  by  other  researchers. 

The  research  was  performed  in  four  coordinated  steps:  1)  solution  of  microscale  problems;  2) 
development  of  homogenization  schemes  for  multiscale  analyses;  3)  construction  of  local  and 
global  damage  evolution  algorithms;  and  4)  solution  of  model  demonstration  problems.  As  a 
part  of  the  effort,  a  new  simplified  explicit  finite  element  code  is  being  developed  by  the 
principal  investigator.  This  code  is  identical  in  general  formulation  to  DYNA3D,  but  is  designed 
to  simplify  the  code  development  process,  and  to  streamline  the  implementation  of  the  newly 
developed  multiscale  subroutines  directly  into  DYNA3D.  All  of  these  results  are  included  in  the 
Ph.D.  dissertation  of  Cha  Searcy,  a  graduate  student  who  has  worked  on  the  project  with  Dr. 
David  Allen,  and  they  are  not  reported  here. 

The  results  of  a  damage  constitutive  model  for  polymer-matrix  composite  materials  that  was 
developed  and  implemented  into  a  commercially  available  finite  element  package  are  discussed 
in  this  report. 


INTRODUCTION 


Continuous  Damage  Mechanics  (CDM)  is  a  more  rational  approach.  PMCs  experience  damage 
and  unrecoverable  deformations  similar  to  plasticity.  CDM  can  be  applied  to  the  constituents 
(fiber,  matrix,  interphase)  at  the  micro-scale  [1-3].  In  principle,  these  models  should  need  a  small 
number  of  material  parameters  because  each  phase  is  isotropic  [3].  However,  the  constituents 
cannot  be  tested  in  isolation,  so  they  must  be  tested  as  part  of  a  lamina.  It  is  difficult  to  design 
material  tests  at  the  mesoscale  (lamina)  level  that  reveal  the  material  parameters  in  the 
constituents.  One  such  micro  scale  level  model  is  presented  in  [2].  The  micro  scale  level  model 
in  [3]  includes  a  detailed  identification  of  material  parameters  in  terms  of  lamina  level  tests. 
From  a  computational  point  of  view,  micro  scale-level  models  are  expensive,  because  they 
require  an  excellent  micro  mechanics  model  in  order  to  assemble  the  contributions  of  the 
constituents  into  the  lamina  behavior.  Furthermore,  accurate  micro  mechanics  is  needed  to 
decompose  the  meso  scale  level  stress  and  strain  among  the  constituents.  Ref.  3  uses  the 
periodic  microstructure  model  (PMM)  from  [4],  which  is  very  accurate. 

A  simpler  approach  is  to  work  directly  at  the  mesoscale  modeling  the  behavior  of  a 
homogeneous,  orthotropic  lamina.  The  material  is  assumed  to  be  orthotropic  and  not 
transversely  isotropic  because  even  isotropic  materials  become  orthotropic  after  damage  [5]. 
Therefore,  use  of  scalar  damage  is  not  appropriate,  although  it  has  been  done  [6].  The  simplest 
yet  realistic  representation  of  damage  is  a  second  order  diagonal  tensor  with  principal  directions 
aligned  with  the  material  directions  [7].  This  works  well  for  PMC  because  the  fiber  direction 
dominates  the  response;  all  micro  cracks  are  thus  oriented  in  three  perpendicular  planes 
coinciding  with  the  material  directions.  Any  more  generality  yields  models  too  complex  to  be 
practical. 
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A  simple  meso  scale  level  model  of  damage  and  plasticity  is  presented  in  [8],  It  works  only  for 
plane  stress  without  inter-laminar  damage.  Its  main  drawback  is  its  cumbersome  identification. 
In  order  to  identify  the  material  parameters,  one  has  to  test  a  [+/-67]s  laminate,  which  is  not 
standard.  Furthennore,  the  model  does  not  necessarily  predict  lamina  failure  in  agreement  with 
any  failure  criteria  such  as  Tsai-Wu. 

In  order  to  address  these  shortcomings,  a  new  damage  constitutive  equation  is  presented  in  [7]. 
Then,  Ref.  [9]  extended  the  damage  model  of  [7]  to  laminates,  perfonned  identification  of 
several  PMC  materials  from  published  data,  and  validated  the  model  with  additional  data  not 
used  in  the  identification  process.  Later,  Ref.  [10]  added  plasticity  to  the  model,  identified  new 
materials  and  further  validated  the  plasticity  and  the  damage  portions  of  the  model  with  new 
data.  The  analysis  is  perfonned  for  a  lamina  in  [7],  and  for  laminates  in  [9-10]  using  CLT. 
Then,  Lonetti  et  al.  [11]  added  inter-laminar  effects.  A  finite  element  formulation  of  the  3D 
damage  model  of  [11]  is  presented  and  implemented  here  as  a  user  material  model  within 
ANSYS  in  order  to  solve  more  complex  laminated  problems  with  general  boundary  conditions 
and  including  the  effect  of  geometrical  non  linearity. 

FINITE  ELEMENT  FORMULATION 

A  displacement-based  finite  element  formulation  is  used.  The  body  of  a  laminate  is  represented 
by  a  series  of  layers  with  different  thickness  and  orientations.  Two  reference  systems,  global  (x,, 
i=1..3)  and  material  (<?,,  i=1..3)  are  used.  The  local  reference  axis  ei  is  aligned  with  the  fibers 
direction,  e3  points  through  the  thickness  of  the  laminate,  and  C2  lays  on  the  mid-surface  of  the 
layer  and  is  perpendicular  to  ei  and  e3  (e2=  e3  x  ei). 


Considering  a  body  B,  in  which  the  internal  stresses  a,  the  distributed  load  q,  and  the 
concentrated  loads/ constitute  an  equilibrated  system,  applying  an  arbitrary  virtual  displacement 
pattern  Sw*  compatible  with  the  internal  strain  8s*,  the  principle  of  virtual  displacement  can  be 
written  as 


(1) 


where  nf  is  the  number  of  the  concentrated  loads.  Introducing  the  shape  functions  N  and  nodal 
displacement  d  and  using  the  finite  element  discretization  process,  Eq.(l)  becomes 

f  \ 


Sd 7 


J  BTadV  - 1 NT qdS  -  f 


=  0, 


(2) 


W  s  J 

where  S  is  the  area  in  the  X1-X2  plane,  and  V  =  [S  x  t\,  with  t  total  thickness  of  the  laminate.  To 
predict  the  stiffness  reduction,  an  incremental  step-by-step  analysis  is  adopted.  For  a  load 
increment,  Eq.(2)  becomes 


|  Bl  AadV  - 1  NT  AqdS  -  2>/  =  0  (3) 

V  S  <=1 


where  A a  represents  the  stress  increment.  The  damaged  stiffness  E  at  the  current  increment  can 
be  detennined  directly  from  the  constitutive  equations  (Eq.  42-43)  as 


.  da  da  ■ 

a  = - s-\ - D  =  Es 

ds  dD 
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Substituting  Eq.  (  4)  into  Eq.(3),  we  get 


where  the  stiffness  matrix  and  the  load  vector  are 

[*J  =  J  [BTEepdB]iv 

V 

[AT7]  =  |  NT  AqdS  +  A/ 

The  incremental  form  of  the  governing  equation  can  be  written  as 

[<+<][  Au]  =  AF  (7) 

where  A u  is  the  incremental  displacement,  K^L  is  the  non  linear  material  tangent  stiffness 
matrix,  K%  is  non  linear  geometrical  contribution,  and  A F  is  out  balance  force. 

Damage  is  monitored  at  the  Gauss  integration  points.  The  following  algorithm  is  used  to 
integrate  numerically  the  rate  equations.  First,  compute  the  strain  and  stress  increment  As,  Agin 
the  local  coordinate  system  for  each  lamina  at  each  gauss  point 

|>4=7;[A4  ;  [A<t]£  =  Tk  [Ao-]g  (8) 


(5) 


(6) 


where  the  Tk  is  a  coordinate  matrix  transformation  [1].  Subscripts  L,  G,  indicate  local  and  global 
coordinates,  respectively.  An  elastic  predictor  and  inelastic  corrector  scheme  is  used  to 
determine  the  effect  of  a  small  strain  increment  As  .  In  this  way  the  initial  increment  is  purely 
elastic.  Damage  is  evaluated  in  order  to  check  if  the  inelastic  effects  grow.  There  are  two 
possible  cases,  elastic  behavior  ( gd  <  0  )  or  damage  evolution  ( gd  >  0  ).  The  evolution  of 
damage  variables  is  subjected  to  the  return-mapping  algorithm.  In  this  way,  the  damage  domain 
is  linearized  to  the  first  order  as 


(y-£')+td  b-A) 


i+ 1 


dy 


/+! 


(9) 


where  the  subscript  (7+1)  indicates  load  step,  while  superscript  (k)  represent  the  iteration  number. 
The  thennodynamic  forces  (Y,y)  and  the  effective  stress  cr ,  can  be  expressed  by  using  the 
constitutive  and  evolution  equations  in  tenns  of  the  kinematic  quantities  (D, s)  and  the  damage 
multiplier  as 
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(4:'-A)=§‘  (W-DLb f  ‘  (4t-4,) 

i+l  ^  i+l 

/  i+l  k  \  dy  A  jd 

\ri+i  -rM)  =  —  -aa 

5^  <+1 

(5!r-ss)=M-'(<c,>-^,>)  no) 

Kr -CT,™) =[s(i>)(+r  -4,)+^‘  K1  -At,) 

i+l 

(Atr-At,)=AA|^‘  (.-r-o 

i+l 

In  the  general  case  when  damage  grows,  Eqs  (10)  leads  to  a  linear  system  that  can  be  solved  for 
the  damage  multiplier  A/l‘/ .  The  last  step  is  to  check  if  the  damage  state  variables  D,  have 
reached  the  critical  values  D] R  .  When  a  damage  component  exceeds  the  critical  value,  it  signals 

the  appearance  of  a  macro-crack.  Therefore,  the  damage  component  is  set  to  nearly  one  ( Dj=l ) 
to  enforce  total  reduction  of  stiffness  at  the  Gauss  point. 

USER  DEFINED  MATERIAL  MODEL 

In  order  to  incorporate  material  non-linearity,  a  user-defined  subroutine  is  written  in  Fortran  and 
linked  with  ANSYS.  In  the  process  of  linking,  a  customized  ANSYS  executable  file  is  obtained 
which  is  used  for  the  analysis.  The  procedure  for  getting  the  executable  file  [12]  is  as  follows:  A 
new  directory  is  created  in  the  drive  where  ANSYS  is  installed.  The  following  files  are  then 
copied  to  the  new  directory  from  the  \custom\user  sub  directory  in  ANSYS:  Anscust.bat, 
Makefile,  Ansysex.def.  The  Fortran  files  are  then  compiled  and  linked  with  ANSYS  by  running 
the  Anscust.bat  file.  The  procedure  will  load  object  files  and  library  files  after  the  compilation 
and  a  new  executable  ANSYS  file  will  be  created  and  will  reside  in  the  new  directory. 

Continuous  damage  mechanics  coupled  with  thermodynamics  is  used  in  order  to  predict  the 
inelastic  behavior  of  a  material.  In  this  case,  a  set  of  internal  variables  is  used  to  describe  the 
damage  behavior.  The  damage  model  used  for  modeling  the  user  material  is  based  on  available 
data  [7,9-1 1].  The  simplicity  of  the  model  lies  in  the  fact  that  only  a  few  parameters  are  required 
to  model  non-linearity  and  they  can  be  obtained  from  standard  tests.  The  model  predicts  the  non¬ 
linear  effects  as  a  reduction  of  stiffness  and  increments  of  damage.  It  requires  20  user  constants, 
including  five  initial  transversely  isotropic  properties,  three  state  variables  di,  do  and  d3,  also 
called  damage  variables,  and  12  material  parameters.  The  damage  variables  are:  the  damage  in 
fiber  direction  d\,  the  damage  in  transverse  direction  th,  and  the  damage  in  the  direction 
perpendicular  to  the  fiber  d3. 

The  12  material  parameters  that  are  required  for  the  damage  model  are  as  follows: 
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a.  Six  internal  parameters  J\  i,  J22,  A,  A,  Ih,  A  that  are  univocally  defined  in  tenns  of 
experimental  material  properties. 

b.  Three  critical  damage  values  in  tension  and  compression  along  the  fiber  and  transverse 
directions:  Atcr,  Accr,  Atcr- 

c.  Three  damage  threshold  and  damage  evolution  parameters:  jc\,  and  Co. 

Internal  Material  Constants 

The  internal  material  constants  (J\  \ ,  J22,  J33,  A,  A,  A3)  are  computed  a  priori  and  do  not  change 
during  the  course  of  the  analysis.  They  are  defined  in  term  of  the  following  experimental 
properties. 

a.  Critical  damage  values  in  tension  (Atcr,  Atcr)  and  compression  (Accr) 

b.  The  stiffness  values  (A,  A,  Gn,  G23,  V12) 

c.  The  strength  values  in  tension  (At,  At),  compression  (Ac)  and  shear  (A,  A,  A) 

d.  Damaged  shear  modulus  at  failure  ( G*n ,  G*n ,  G2*3 ) 

The  internal  constants  are  defined  by  fourth  and  second  order  tensors  J  and  H.  They  appear  in 
the  formulation  of  the  damage  surface  g 

g = (Y,  JpkYhl  f-  +  (|h„  y.jI)1'2  -  r(5)  -  r,  an 


in  thermodynamic  force  space  Y 

i  S(E„rle„e„) 

^  2  DDy  ( 

In  stress  space,  Eq.  (11)  represents  the  Tsai-Wu  surface  at  failure  with 
Y  =  Thennodynamic  force  tensor  dual  to  the  damage  tensor  D 
J,  H  =  Internal  material  constants 
y  (5)  =  Damage  evolution  variable 

Yo  =  Damage  threshold  representing  the  initial  size  of  the  damage  surface 
No  damage  occurs  until  the  thermodynamic  force  Y  reaches  the  damage  surface.  For  undamaged 
material  y  =  0  and  g  has  the  shape  of  the  Tsai-Wu  surface. 

At  failure,  y  +  yo  =  1  and  the  shape  and  size  of  g  matches  the  Tsai-Wu  surface,  where  y 
represents  the  value  of  y  at  failure.  Comparing  the  g-surface  with  the  Tsai-Wu  surface  at  failure, 
we  arrive  at  a  linear  system  of  equations  for  detennination  of  the  J  and  H  tensors  as  follows. 


Calculation  of  Jn,  Hi 

When  the  composite  lamina  is  subjected  to  uniaxial  load  in  the  longitudinal  direction,  all  stress 
components  other  than  oi  are  zero.  Therefore  (Eq  .11)  reduces  to 


=  /r^ 

^ 11  a6 


N^r0'  -  (t+a) 


with  the  integrity  tensor  given  by 
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Q  =  yjl-D 


(14) 


where 


C  =  Compliance  tensor  of  the  undamaged  material 
Q  =  Integrity  tensor 
D  =  Damage  tensor 

At  failure,  the  damage  variables  reach  the  critical  values  (Du,  D ic).  If  Fit  represents  the  tensile 
strength  and  Flc  the  compressive  strength,  (Eq.  13)  becomes 


H.|§tFm  =  (t‘+To) 

W.lSrF.c  =  (t  +y0) 

iiic 

and  Qlc  =  1-Dlc 


(15) 

(16) 
(17) 


The  Tsai-Wu  criterion  for  uniaxial  loading  in  the  fiber  direction  is  given  by 

1  and  f.F,  +  f , ,  F,2.  =  1 


f  F  +f  Fz 

iirit  unit 


(18) 


Hence  the  right  hand  side  of  (Eq.15)  and  (Eq.16)  should  equal  1  so  that  the  damage  surface 
matches  with  the  Tsai-Wu  surface  at  failure.  The  critical  damage  values  are  obtained  from 
statistical  methods  [7].  Then  the  two  equations  are  solved  simultaneously  and  the  values  of  Jn 
and  Hi  are  obtained. 


Calculation  of  J22,  H2 

When  the  composite  lamina  is  subjected  to  transverse  uniaxial  loading,  all  the  stress  components 
other  than  02  are  zero.  At  failure,  (Eq  .11)  reduces  to 


(19) 


D2t  1  D2t 


(20) 


When  the  lamina  is  subjected  to  in-plane  shear  loading,  all  the  stress  components  other  than  05 
are  zero.  At  failure,  (Eq  .11)  becomes 
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'  H,  H2 

2C66 

^  ^2s 

^2s 

(y*  +Yo)  =  1 


(21) 


Since  the  shear  response  of  the  lamina  does  not  depend  on  the  sign  of  the  shear  stress,  the 

coefficient  of  the  linear  term  in  (Eq.21)  should  be  zero.  Therefore  we  get 

*2 


H2  = 


Q.I 

— y-Hj  =  ■ 


r  IT  .  r  _  2s 
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Now  ksl2  can  be  approximated  as  the  ratio  of  damaged  shear  modulus  to  the  undamaged  shear 
modulus  [7,9-1 1] 
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Also,  it  has  been  experimentally  observed  that  most  of  the  shear  damage  is  in  the  form  of 
longitudinal  cracks  rather  than  the  transverse  cracks  [7]  and  hence  D2S>Dis  or  from  (Eq.  12)  we 
obtain  the  restriction 

0  <  rsu  <  1  (25) 

Substituting  Eq.  (24)  in  Eq.  (23)  and  solving  Eq.  (19,  22,  23)  we  obtain  the  values  of  J22  and  H2. 

Calculation  of  J33  and  H3 

In  this  case  the  inter-laminar  stresses  are  taken  into  consideration.  The  formulation  of  equations 
is  similar  as  that  of  the  in-plane  case.  When  the  lamina  is  subjected  to  inter-laminar  stresses,  at 
failure,  (Eq.  1 1)  reduces  to 
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Since  the  shear  response  does  not  depend  on  the  sign  of  the  shear  stress,  the  coefficients  of  the 
linear  term  must  be  zero.  Therefore  we  get, 

H,=-§fH,= ;  r.„=S-  (28) 
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Also,  it  has  been  experimentally  observed  that  rsl3  should  be  less  than  1  [11].  Similar  to  ksl2  in 
(Eq.  24)  we  have  ksl3  and  ks23  given  by 

k.,!=n^L=§i  (30) 
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Therefore  (Eq.  26)  and  (Eq.  27)  reduce  to 
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F4  =  (y*  +Yo)  =  1 


Solving  (Eq.  32)  and  (Eq.  33)  we  obtain  the  values  of  J33,  H3.  The  next  step  is  to  determine  the 
evolution  or  hardening  parameters  ci,  C2  and  damage  threshold  yo. 
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Flow  and  Hardening  Rules 

A  non-associated  flow  rule  is  used  for  the  damage  model  [7,9-1 1].  The  flow  potential  surface  is 
given  by  the  following  equation 

/  =  (ys  JiJllYhk)”-y(5)-To  (34) 


The  damage  and  flow  surface  expand  as  a  function  of  evolution  variable  y  according  to  the 
evolution  law 
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where  7t(S)  is  the  dissipation  energy  and  Ci,  C2  are  material  constants,  which  have  to  be 
determined  from  experimental  data.  Since  the  dissipation  energy  should  be  convex  [7],  its  second 
derivative  should  be  positive 
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In  that  case  the  signs  of  c i,  C2  should  be  different.  The  flow  rule  for  damage  and  hardening  is 
given  by 

dD  =  dX  —  ;  dS=dX—  (37) 

dY  dy 


where  dX  is  the  damage  multiplier  whose  value  can  be  detennined  from  the  consistency 
condition  (g  =  0,  dg  =  0) 
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Substituting  in  (Eq.  37)  we  obtain 
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The  incremental  stress-  strain  relations  for  damage  evolution  is  given  by 
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(40) 
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where  s  is  the  total  strain.  Substituting  dD  from  (Eq.  39)  into  (Eq.  40)  we  get, 
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Here  Emni  represents  the  tangent  stiffness  due  to  the  material  non-linearity.  When  geometrical 
non-linearity  is  considered  (Eq.  41)  is  written  as 

d(J  =  (E,nnl  +  Egnl)d£  (43) 


Hardening  Parameters  and  Damage  Threshold 

The  hardening  parameters  (ci,  C2)  control  the  damage  evolution  and  the  damage  threshold  (yo) 
represent  the  initial  size  of  the  damage  surface.  Since  the  material  behavior  is  highly  non-linear 
for  a  composite  lamina  in  the  in-plane  shear  mode,  the  damage  is  assumed  to  be  notable  in  this 
case  [7].  Therefore,  ci,  Co  and  yo  are  adjusted  to  predict  the  shear  response  of  the  lamina 
subjected  to  pure  shear  conditions  using  Finite  Element  Analysis. 

A  single  lamina  is  modeled  in  ANSYS  and  subjected  to  pure  shear.  The  material  properties  and 
parameters  are  input  in  the  ANSYS  material  model  definition.  The  rear  nodes  are  clamped,  and 
the  nodes  in  the  side  faces  are  free  to  move  only  in  the  y-direction.  The  nodes  in  the  front  face  of 
the  lamina  are  given  a  unifonn  displacement  v  in  the  y-direction  to  simulate  the  pure  shear 
condition. 

The  non-linear  analysis  is  run  with  the  user  material  model,  which  includes  ci,  C2  and  yo.  During 
the  post  processing  stage,  the  sum  of  the  reaction  forces  Fxy  and  the  displacement  v  in  the  front 
face  of  the  lamina  are  recorded  for  each  sub  step.  The  average  shear  stress  is  calculated  by 
dividing  Fxy  by  the  shear  area  and  the  shear  strain  is  calculated  from  the  displacement  of  nodes  in 
the  front  face  of  the  lamina.  The  shear  stress-strain  from  the  analysis  is  plotted  and  compared  to 
the  experimental  shear  response  (Fig.  1).  If  the  curves  don’t  match,  then  the  ci,  C2  and  yo  values 
are  adjusted  and  the  procedure  is  repeated  until  the  shear  stress-strain  matches  the  experimental 
shear  response. 


11 


GEOMETRIC  AND  FINITE  ELEMENT  MODELING 


Finite  element  models  of  the  [0/90]s  and  [90/0]s  lay-up  beams  with  different  aspect  ratios  (length 
to  height  ratio)  are  developed.  Four  point  bending  tests  are  used  in-order  to  study  the  effect  of 
aspect  ratios  on  the  damage  of  the  beams.  Both  material  and  geometric  non-linearity  are 
considered  for  the  analysis  [14].  Contributions  due  to  in-plane  and  total  damage  (both  in-plane 
and  inter-laminar  damage)  are  analyzed  for  each  aspect  ratio  of  the  beam. 

Laminated,  rectangular  beams  with  different  aspect  ratios  are  modeled  in  ANSYS  using  solid 
elements.  Modeling  with  shell  elements  is  not  possible  within  ANSYS  because  the  code  does  not 
supply  inter  laminar  stresses  to  the  user  subroutine.  Other  codes  such  as  ABAQUS  and  LUSAS 
were  investigated  and  the  same  limitation  was  found.  Due  to  symmetry,  only  a  quarter  of  the 
beam  is  modeled.  A  rectangular  prism  is  modeled  based  on  the  aspect  ratio  of  a  beam  and  then 
partitioned  into  4  layers  representing  the  lay-up  configuration  with  different  orientations.  Three 
aspect  ratios  are  considered  (a/h  =4,  10,  and  20).  Then,  the  models  are  meshed  using  3-D  20- 
noded  structural  solid  elements.  The  material  properties  vary  with  the  orientation  of  the  fibers  in 
each  layer.  Therefore  the  material  orientation  option  is  used  in-order  to  orient  the  elements,  in 
each  layer,  according  to  the  lay-up  angle.  Orthotropic  material  properties  based  on  available  data 
[7,9-11]  are  assigned  to  each  layer.  Glass/Epoxy  and  Carbon/Epoxy  materials  are  considered  for 
the  analysis.  The  material  properties  are  shown  in  Table  1. 


Since,  only  a  quarter  of  the  beam  is  considered,  symmetric  boundary  conditions  are  given  to  the 
nodes  at  the  mid  span  on  the  plane  YZ  and  also  on  the  plane  of  bending  XZ  (see  Fig.  8.3  in  [1]). 
The  nodes  in  the  mid-plane  at  the  end  of  the  beam  are  constrained  in  the  z-direction  to  simulate 
the  end  support.  The  nodes  at  a  distance  of  l/3rd  of  the  length  of  the  beam  from  the  ends,  which 
also  belong  to  the  mid  surface  are  given  a  specified  displacement  z-direction  to  simulate  the 
loading.  These  boundary  conditions  simulate  a  four  point  bending  test.  The  material  parameters 
used  for  the  analysis  are  shown  in  Table  2.  Each  material  is  analyzed  in  the  [0/90]s  and  [90/0]s 
lay-ups.  The  material  properties  and  parameters  (with  adjusted  ci,  C2  and  yo)  are  input  in  the 
ANSYS  material  model  definition.  The  analysis  is  performed  considering  (a)  In-plane  damage 
only,  and  (b)  Total  damage  (in-plane  and  inter- laminar  damage). 


In-plane  Damage  Only 

In  this  case,  the  material  parameters  J33  and  H3  are  set  to  zero  so  that  the  inter-laminar  damage  is 
turned  off.  The  non-linear  analysis  is  performed  on  the  beam  with  the  optimum  number  of  sub 
steps.  The  damage  is  monitored  at  the  integration  points  of  the  elements.  When  the  damage  state 
variables  di  reach  the  critical  values  D’r ,  the  damage  component  is  set  to  1  in  order  to  enforce 

total  stiffness  reduction  at  the  Gauss  point.  During  post  processing,  the  sum  of  reaction  forces  R~ 
at  the  restrained  nodes  and  the  deflection  w  of  the  end  nodes  located  in  the  mid-plane  are  noted. 
Then  the  bending  stress  S  is  calculated  in  tenns  of  the  reaction  R-  output  from  the  analysis  as 


_  Me  _  2aR„ 


(44) 


where  a,  b,  h  are  the  length,  width  and  height  of  the  beam  and  R-  is  the  reaction  force  measured. 
Then  the  bending  stress  is  normalized  with  respect  to  the  transverse  tensile  strength  of  the 
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material  F2t.  Similarly,  the  deflection  w  is  normalized  with  respect  to  the  length  of  the  beam  a.  A 
graph  is  plotted  between  and  w/a  of  the  beam.  This  procedure  is  repeated  for  all  types  of 
beams  (lay-ups  and  materials).  In  case  the  deflection  of  the  beam  becomes  large  (more  than  30% 
of  the  beam  length),  geometric  non-linearity  is  considered  in  addition  to  the  material  non¬ 
linearity. 

Total  Damage 

In  this  case,  the  inter-laminar  damage  is  turned  on  when  actual  values  of  J33  and  H3  are  input. 
The  rest  of  the  procedure  is  the  same  as  the  previous  case.  The  evolution  of  the  damage  variables 
da,  d3  (state  variables)  is  studied.  The  maximum  value  of  the  damage  variables  is  noted  at  the 
integration  points  for  each  sub-step  and  a  graph  is  plotted  between  d2,  d3  and  w/a.  The  graph 
gives  an  indication  of  the  damage  growth  in  the  beam. 

RESULTS  AND  DISCUSSION 

The  damage  threshold  yo  and  hardening  parameters  ci,  C2,  are  adjusted  to  obtain  a  good 
correlation  between  the  model  predictions  and  experimental  shear  stress-strain  data,  as  illustrated 
in  Fig.  1 .  The  results  obtained  for  the  Glass/Epoxy  composite  beam  with  different  aspect  ratios 
(a/h  =  4,  10,  20)  and  lay-up  angles  ([0/90]s  and  [90/0]s)  are  shown  in  Figs.  2-10.  With  reference 
to  Fig.  2,  the  straight  line  indicates  that  there  is  no  damage  considered  for  the  analysis.  This 
condition  is  achieved  by  assigning  a  large  value  for  the  damage  threshold  (yo=lE20),  which 
means  that  the  initial  damage  surface  is  very  large  and  hence  no  damage  occurs.  The  curve 
below  the  linear  case  is  the  one  obtained  by  considering  the  in-plane  damage  only.  The  main 
mode  of  damage  is  due  to  d2,  as  shown  in  Fig.  4.  When  d2  reaches  the  maximum  value,  the 
solution  does  not  converge  and  the  program  stops.  This  indicates  the  emergence  of  a  macro 
crack.  The  maximum  damage  value  d2  is  0.59,  which  exceeds  the  actual  critical  damage  value 
D2cr  (Table  2).  This  is  because  the  damage  value  0.59  occurs  only  in  a  small  region,  as  it  also 
happens  for  d3  (see  Fig.  3.)  When  D2cr  is  considered  as  the  limiting  damage  value,  the  analysis 
stops  when  d2  reaches  D2cr,  as  shown  in  Fig.  2. 

The  last  curve  indicates  the  results  for  which  the  total  damage  is  considered  (in-plane  and  the 
inter-laminar  damage).  In  this  case  d3  dominates  the  damage  behavior.  This  is  because  the 
Glass/Epoxy  is  weak  in  the  inter-laminar  direction.  This  is  further  shown  in  Fig.  4  where  it  can 
be  seen  that  d3  is  much  higher  than  d2.  The  values  plotted  are  the  maximum  values  of  the  damage 
variables  obtained  from  the  integration  points  in  each  sub-step.When  D3Cr  is  considered  as  the 
limiting  damage  value,  the  analysis  stops  when  d3  equals  D3cr,  as  indicated  in  Fig.  2.  The 
maximum  damage  value  (ds)  is  0.535,  which  exceeds  the  actual  critical  damage  value  D3Cr  (Table 
2).  The  reason  being  is  that  d3  equals  0.535  at  a  small  region,  as  shown  in  Fig.  3. 

It  is  worth  noting  that  for  the  [0/90]s  laminate,  the  maximum  damage  in  both  in-plane  (d2)  and 
total  damage  (d3)  condition  occur  in  the  90°  laminate.  This  is  because  the  0°  laminate  is  subjected 
to  tension  at  the  bottom  and  compression  at  the  top,  where  damage  is  very  small  because  the 
fibers  are  strong  in  the  0°  direction.  At  a/h=4  (Fig.  2-4),  inter-laminar  shear  is  high  and  it  causes 
significant  out  of  plane  damage  d3.  If  the  model  does  not  consider  it  (d3=0),  the  in-plane  damage 
is  underestimated  (Fig.  4).  When  the  aspect  ratio  is  20,  the  damage  d2  becomes  significant  in  the 
total  damage  condition,  as  shown  in  Fig.  5-6.  The  maximum  values  of  d2  (0.38)  and  d3  (0.46) 
become  closer,  which  means  that  d2  begins  to  dominate  the  damage  in  the  beam  for  higher  aspect 
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ratios.  This  is  because  inter-laminar  shear  is  less  pronounced  at  high  aspect  ratios.  The  effect  of 
damage  on  deflections  (Fig.  5)  for  a/h=20  is  small  because  the  outer  layers  dominate  the  stiffness 
of  the  beam.  Still,  the  inter-laminar  shear  strength  F4  of  Glass/Epoxy  is  weak,  resulting  in 
significant  inter-laminar  damage  d3  (Fig.  6). 

Similar  results  can  be  observed  in  the  [90/0]s  lay-up  beams  (Fig.  7-10).  As  a/h  increases,  the 
mode  of  failure  changes  from  d3  to  d2  (Table  3).  In  the  case  of  the  [90/0]s  laminates,  even  for  the 
aspect  ratio  of  10,  d2  starts  to  dominate  the  damage  in  the  total  damage  condition.  This  is  due  to 
the  fact  that  the  90°  laminate  is  subjected  to  a  high  stress  at  top  and  bottom  of  the  beam  in  the 
transverse  direction  where  the  material  is  weak.  This  is  the  reason  why  the  maximum  value  of  d2 
occurs  in  the  90°  laminate  and  the  maximum  value  of  d3  occurs  in  the  0°  laminate.  When  the 
aspect  ratio  is  20,  very  high  deflections  in  the  beam  are  observed.  Hence,  geometrical  non-linear 
effects  must  be  considered  (Fig.  9).  The  mode  of  failure  switches  to  transverse  damage  d2  at 
a/h=20  (Table  3  and  Fig.  10). 

The  inter-laminar  damage  d3  dominates  the  damage  in  all  cases  for  Carbon/Epoxy  composite 
beam  with  different  aspect  ratios  (a/h  =  4,  10,  20)  and  lay-up  angles  ([0/90]s  and  [90/0]s).  In  case 
of  the  [0/90]s  lay-up  beams  (Fig.  1 1-12),  d2  is  significant  if  d3  is  not  modeled  (labeled  as  d3=0  in 
the  figures),  whereas  when  inter-laminar  effects  are  considered,  d3  controls  the  damage  in  the 
beam.  This  is  due  to  the  fact  that  the  material  is  strong  in  the  transverse  direction  and  withstands 
damage  due  to  d3  even  when  the  aspect  ratio  increases  and  the  in-plane  damage  is  insignificant. 

For  Carbon/Epoxy  [0/90]s,  the  dominant  mode  of  failure  remains  d3  even  at  a/h=20  (10-11).  The 
effect  of  damage  on  deflections  is  negligible  because  the  fibers  on  the  outer  layers  dominate  the 
bending  stiffness.  In  case  of  the  [90/0]s  lay-up  beams,  similar  results  as  the  previous  case  is 
obtained  except  that  geometric  non  linear  effects  are  also  considered  here  due  to  the  large 
deflection  of  the  beam.  Even  switching  LSS  to  [90/0]s  the  mode  of  failure  remains  inter-laminar 
damage  d3  for  Carbon/Epoxy  with  4<a/h<20  (Table  4).  Significant  transverse  damage  d2  is 
accumulated  but  it  is  not  the  critical  factor.  Neglecting  d3  leads  to  gross  overestimation  of  the 
failure  load. 


CONCLUSIONS 

A  finite  element  formulation  of  continuum  damage  model  that  includes  inter  laminar  effects  is 
successfully  developed  and  implemented  as  a  user  subroutine  in  a  commercial  finite  element 
analysis  code.  The  internal  damage  parameters  are  identified  from  available  experimental  data. 
Furthermore,  the  damage  threshold  and  evolution  parameters  are  identifies  from  available, 
experimental  shear  stress  strain  data.  The  finite  element  formulation  and  identifies  material 
parameters  are  then  used  to  investigate  the  influence  of  material  properties,  aspect  ratio,  laminate 
stacking  sequence,  and  geometric  non  linearity  on  the  damage  evolution,  deflections,  and  failure 
of  laminated  beams  under  four  point  loading.  Furthennore,  the  option  of  neglecting  inter  laminar 
damage  is  shown  to  lead  to  misleading  predictions  of  in  plane  damage  evolution  and  poor 
predictions  of  failure  load.  In  case  of  Glass/Epoxy  composite,  inter-laminar  damage  is  significant 
for  all  [0/90]s  and  [90/0]s  beams  with  low  aspect  ratios.  As  the  aspect  ratio  increases,  the  in-plane 
damage  becomes  dominant.  In  the  case  of  Carbon/Epoxy  composite,  the  inter-laminar  damage  is 
dominant  for  all  aspect  ratios  of  the  beam. 
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SIGNIFICANCE  -  ARMY  VALUE 


Weapons  technology  is  shifting  increasingly  towards  minimizing  mass  of  defensive  weapons  as  a 
primary  design  constraint.  This  is  readily  apparent  in  both  tank  armor  and  single  soldier  armor, 
not  to  mention  numerous  other  defense  applications.  As  a  part  of  this  new  design  philosophy,  it 
is  clear  that  significant  decreases  in  armor  weight  can  be  achieved  by  utilizing  multi-layered  and 
multi-material  structures  to  defeat  offensively  fired  weapons. 
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Table  1.  Material  properties  of  Glass/Epoxy  and  Carbon/Epoxy 


Property 

Glass/Epoxy 
(Vf=  52%) 

Carbon/Epoxy 
(Vf=  62%) 

Ei  (MPa) 

36200 

131300 

E2  (MPa) 

8100 

10800 

G 12  (MPa) 

3200 

5200 

Gi3  (MPa) 

1500 

7100 

Vl2 

0.23 

0.29 

Fit  (MPa) 

932 

2100 

Fic  (MPa) 

370 

1080 

F2,  (MPa) 

42 

80 

F4  (MPa) 

35 

72 

F5  (MPa) 

50 

70 

F6  (MPa) 

50 

70 

G12  damaged  (GPa) 

2700 

2258 

G13  damaged  (GPa) 

1299 

3083 

G23  damaged  (GPa) 

1299 

3083 

Table  2.  Material  Parameters  for  the  damage  model 


Property 

Glass/Epoxy 

Carbon/Epoxy 

Jll 

0.53  le-2 

0.1897e-2 

J22 

0.226 

0.287e-l 

J33 

0.2009 

0.534e-l 

HI 

0.681e-l 

0.2571e-l 

H2 

-0.1445e-l 

-0.807e-2 

H3 

0.592e-l 

0.8018e-2 

ks12 

0.866 

0.434 

ks13 

0.866 

0.434 

ks23 

0.866 

0.434 

rsi2 

0.212 

0.314 

rsn 

0.8703 

0.311 

rs23 

4.105 

0.992 

Cl 

0.22 

0.045 

c2 

-0.3 

-0.625 

Yo 

0.015 

0.07 

D!cr 

0.1161 

0.1161 

D2cr 

0.5 

0.5 

D3cr 

0.5 

0.5 
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Table  3.  Failure  stress  S/F2t  and  mode  of  failure  for  Glass/Epoxy 


Laminate 

stacking 

sequence 

a/h 

S/F2t 

Mode 

At  Failure 

d3=0 

d#0 

[0/90]g 

4 

1.57 

1.22 

d3 

d2=0.243 

10 

4.78 

3.61 

d3 

d2=0.252 

20 

10.00 

7.76 

d3 

d2=0.375 

[90/0]s 

4 

4.12 

1.01 

d3 

d2=0.204 

10 

5.04 

2.91 

d3 

d2=0.493 

20 

3.45 

1.79 

d2 

d3=0.283 

Table  4.  Failure  stress  S/F2t  and  mode  of  failure  for  Carbon  Epoxy 


Laminate 

stacking 

sequence 

a/h 

S/F2t 

Mode 

At  Failure 

o 

ll 

d3#0 

[0/90]s 

4 

1.61 

0.69 

d3 

d2=0.0058 

10 

5.30 

2.12 

d3 

d2=0.0061 

20 

11.42 

5.68 

d3 

d2=0.0069 

[90/0]s 

4 

5.63 

0.50 

d3 

d2=2.4E-6 

10 

5.53 

1.67 

d3 

d2=0.2271 

20 

3.20 

2.00 

d3 

d2=0.3694 
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Fig.  3  Damage  evolution  CI3  for  total  damage  analysis  of  [0/90]s  Glass/Epoxy,  a/h=4,  at 
w/a=  0.00913  and  normalized  stress  S/F2t=1.219 
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